ISEA Inverse projection (#3047) #4211
Conversation
|
@rouault @MichaelJendryke @ppKrauss @ldesousa would you be able to help test / validate? Thanks! |
jerstlouis
force-pushed
the
master
branch
5 times, most recently
from
e734459
to
f6805d9
Compare
|
This inverse implementation is now working (successful latlon==>ISEA==>latlon roundtrip) now assuming the parameters:
For For other values of For the Will different values of Not sure about supporting the other parameters. I've clarified that our (provisionary) OGC:1534 CRS is equivalent to: which implies a geodetic->geocentric conversion that does not get performed automatically by Not sure whether this should automatically be done by For the ISEA 5x6 space (provisionally called OGC:153456) used to define the ISEA9R 2D Tile Matrix Set (where 10/30 level 0 tiles are valid -- could also define an ISEA4R quad-tree 2DTMS), which is also a convenient space for doing ISEA3H calculations, we can simply add an affine transform: In the 5x6 space, the twelve icosahedron vertices are very conveniently/intuitively at:
and these points are all on the equator:
|
|
For the appveyor MSVC failure, maybe we need to add |
|
The need for the pipeline to do the correct geodetic->geocentric conversion is problematic with e.g., whereas this kindof works now (complained about missing inverse before): except for some missing bits at the extremities:
|
jerstlouis
force-pushed
the
master
branch
4 times, most recently
from
821caf7
to
bb1609a
Compare
|
@rouault This inverse seems to be working quite well. This is the sample GeoJSON file I mentioned in #3047 inverse-projected to EPSG:4326 in QGIS by doing a Set Layer CRS to:
(together with a Natural Earth countries layer in EPSG:4326). Or de-projecting countries of the whole world in ISEA projection:
Here's how a full ISEA GEBCO raster (https://maps.gnosis.earth/ogcapi/collections/gebco/map.tiff?crs=OGC:1534) deprojects in QGIS (just a few minor glitches around the interruptions):
Should I file a separate issue about whether I believe that the correct ISEA projection actually implies that step, and just as mentioned in OSGeo/gdal#10537, it will still not be possible to accurately visualize things since the cc. @ppKrauss @mpadillaruiz @ldesousa @perrypeterson @allixender |
no idea. What does the literature about this projection says about that? But the fact that you define a spherical ellipsoid with +R=6371007.18091875 and you mention the authalic radius makes me think that perhaps you would need to use the +R_A switch instead (https://proj.org/en/9.4/usage/ellipsoids.html#ellipsoid-spherification-parameters instead), which is sometimes used for a number of projection methods |
|
Thanks for triggering the builds @rouault . The compilation warning as errors settings are a bit over the top... Any chance for some help addressing or bypassing those issues, or is there otherwise an easy way to test this locally? cmake / make / make test all work perfectly fine on my Linux system (g++ (GCC) 14.1.1 20240720), but I don't get those warnings as errors on what seems to be the compilers' opinionated idea of proper C++ style (not being able to trigger the CI builds myself doesn't help either). Is there a particular set of GCC flags that will likely results in all builds passing if I manage to compile locally with them?
As far as I know the official ISEA paper is the Snyder1992 paper. But this paper may actually be describing multiple projections as @ldesousa mentioned, and it definitely mentions multiple polyhedrons such as the truncated icosahedron for example, whereas ISEA is really about the icosahedron base solid (Figure 12, though that does not show the standard ISEA orientation but the The term geodetic or geocentric do not appear, and the word latitude appears on page 14 and 20 only, and the letter ϕ. and lat is used and "geographic laititude" is used which would make you think it it is talking about geodetic latitude. The main reason I've started to think this geodetic ==> geocentric step was necessary is because it was needed to correctly integrate data served by PYXIS / uCalgary using their PYX0 format, for things to align properly in our visualization client. We mention this in the FMSDI Phase 3 ER. I am hoping that experts on the topic (@sahrk @franz-benjamin @shatzi @perrypeterson @mpadillaruiz @ldesousa @allixender) could chime in, so that we could clarify whether the geocentric conversion is implied or not. Essentially we need an agreement on what the ISEA projection really is. This is why a new issue might be a better place for this discussion than this pull request or the original missing inverse issue...
I did try the The #define wgs84InvFlattening 298.257223563
#define wgs84Major 6378137.0
#define wgs84Minor (wgs84Major - (wgs84Major / wgs84InvFlattening)) // 6356752.3142451792955399
#define a2 (wgs84Major * wgs84Major)
#define c2 (wgs84Minor * wgs84Minor)
#define eccentricity 0.0818191908426 // sqrt(1.0 - (c2/a2));
#define log1pe_1me 0.1640050079196 // log((1+e)/(1-e)))
#define S (Pi * (2 * a2 + c2/eccentricity * log1pe_1me))
#define earthAuthalicRadius 6371007.18091875 //sqrt(S / (4 * Pi)); // R -- authalic sphere radius for WGS84 [6371007.1809184747 m]This also matches the numbers here -- √(510 065 621.724 km2 / (4 * π)) = The documentation is a bit confusing, compounded by the fact that both
I imagine this ellipsoid defaults to GRS80? So if I want to use an authalic sphere with the same surface area as WGS84, would I say Omitting the However, I still get a small difference between the Possibly related discussion by other confused users. Does any of this have any relation to the extra Thanks! |
Snyder's paper gives projections for multiple polyhedra using the same approach. For convenience we named his icosahedral projection "ISEA".
The ISEA projection is defined on the sphere. On the sphere there is no distinction between geodetic and geocentric latitudes; they are always equal. If you're integrating data defined on an ellipsoid then some sort of conversion will be necessary. Note that we chose the authalic radius for the ISEA3H DGGS so that authalic latitude could be used to define equal area cells on the ellipsoid. I believe any other conversion approach results in cells on the ellipsoid which are not equal area.
I don't think Snyder specified a "standard" orientation for his projection. Like the acronym ISEA, the orientation we chose for ISEA-based DGGS was something we came up with later.
He's referring to the hexagonal faces of a truncated icosahedron, which is one of the polyhedra he defines his projection on. This is unrelated to the hexagonal cells of a DGGS. |
|
@jerstlouis This is a commendable effort, glad that you took upon it. I can't test it right now, but will keep an eye on this PR. The projection we call "Snyder Equal-area" is actually defined for a right triangle on the spherical surface. It can be used to project a sphere into any polyhedron whose faces can be decomposed into right triangles. Snyder himself was primarily concerned with the Dodecahedron (his favourite), the Icosahedron and the truncated Icosahedron, but the same formulation should apply beyond those. Neither John Snyder, nor those who followed immediately on his steps towards DGGS (White, Kimmerling, etc), defined a "standard" or "default" orientation of the polyhedra. Concerning the Icosahedron, those articles from the early 90s debate to some extent the "inconvenience" of an icosahedron positioned with vertices on the poles always having at least one other vertex on land. In conclusion: the orientation of the polyhedron should be a parameter to any implementation of Snyder's Equal-area projection. That could be the position of a vertex plus an azimuth to a proximate (adjacent?) vertex. That's what @sahrk did in dggrid. |
|
@sahrk Thank you very much for the clarifications!
I guess this means that
Assuming that this authalic latitude conversion is not already done by the ISEA projection in PROJ and geogrid, then I suspect that our implementation of the ISEAH grid as well the PYXIS implementation we tested with in the FMSDI might not be the true ISEA3H DGGS, since we seemed to be using the geocentric latitude rather than the authalic latitude. @perrypeterson @shatzi Any insights / evidence of the contrary about that would be very useful. The zones are clearly not perfectly equal area as computed on the WGS84 ellipsoid (though within the 1% margin allowed by Topic 21 v2 "Equal Area") as can be seen here. I am curious about whether geogrid is actually using the geodetic latitude for its DGGS, or converting from/to geocentric or authalic elsewhere in the code? I could not find anywhere in the code where this is done. @franz-benjamin Any insight on that would be very useful. So ideally this OGC:1534 CRS which we want to define for the purpose of facilitating the implementation of ISEA DGGs should probably imply this geodetic ==> authalic latitude conversion before applying the ISEA projection (rather than geodetic ==> geocentric). That currently seems even more difficult to do with PROJ, since that is not one of the conversions already supported by the pipelines. @sahrk Does this sound like a reasonable assessment? In order to achieve interoperability between multiple implementations of ISEA DGGs, we need clarity on how to bring in and take out georeferenced information correctly, including how geodetic latitude is mapped to the ISEA sphere prior to mapping to the icosahedron. Right now it seems like different implementations of ISEA3H are using three different types of latitudes: geodetic (geogrid), geocentric (PYXIS & GNOSIS), and authalic (DGGRID). NOTE: GDGGS03 does not mention authalic latitudes.
Currently, PROJ supports This implementation of the inverse supports both of these orientations (and could easily support additional ones), but does not yet support the azimuth parameter. I'm also just realizing that we may need an additional parameter reflecting this azimuth for describing the orientation of polyhedron-based DGGs in the draft DGGRS definition -- e.g., to support the R. Buckminster Fuller's Dymaxion Orientation. Thanks all! |
@jerstlouis you're welcome!
The ISEA projection as defined by Snyder and implemented in DGGRID is a spherical projection only. So that sounds logical to me.
I got a 503 error from your data server, but one caveat on cell area variance: if you’re calculating the area on the earth’s surface by treating the edges as great circle arcs then you’re not going to get the exact correct area for an ISEA cell, even on the sphere.
If you want an ellipsoidal datum and you want equal area cells then I think that’s your only choice. btw I’m not a fan of “zones” for cells. It minimizes the important (if so far unrealized) use cases where cell indexes indicate points, rather than areal zones.
I’ll give you a fourth. Probably the most widely disseminated and used dataset on an ISEA grid is the ESA’s SMOS soil moisture products. They wanted 15km cell node spacing, so chose an ISEA4H resolution 9. They also needed it on the WGS84 ellipsoid. So they just took the spherical latitudes produced by DGGRID and treated them as WGS84 ellipsoid latitudes. At a resolution of 15km the error is too small to worry about; it may even be below the representational precision. I’ve referred to that faux ellipsoidal projection as “ISEL”.
I don’t think it mentions ellipsoids/spheres at all; that was beyond the scope of the paper. But DGGRID, which was essentially companion software for that paper, has always used the authalic radius. |
Sorry, the server seems to be prone to these since this week's last update, and we're currently struggling to figure out the cause. I restarted it so hopefully you can try it again before it goes down. Trying to fix these ASAP.
We are not doing that -- we're using several points along the edges and de-projecting each of these to calculate this accurately. Here's a sample image of what that looks like for level 0 zones and the same info as a table below (in summary, compared to the mean zone area (Earth surface area / 12 zones), the 4 polar (Earth geographic poles) level 0 zones have 0.2% smaller areas, whereas the 4 equatorial level 0 zones have 0.17% larger areas, while the remaining 4 zones have 0.03% larger areas):
Level 0 Zones
Total surface area: 510,065,624.74 km²
I am really happy you bring this up! We had some recent discussion for DGGS-JSON where I was advocating for a As for the term "zone", itself, the two terms are defined as such in OGC Topic 21:
So in both cases, they seem to imply this areal meaning.
Wouldn't that be the same as the first geodetic latitude case I mention that I believe geogrid may be using?
Using a sphere with a radius of 6 371 007.1809 m corresonding to √(510 065 621.724 km2 / (4 * π)) (the surface area of the WGS84 ellipsoid) is one thing, which we're doing in both GNOSIS and PYXIS, and which can be specified to PROJ with But using a preliminary conversion of geodetic ==> authalic latitude is a completely different thing, isn't it? i.e., when integrating data referenced to the WGS84 ellipsoid, applying this step to convert geodetic latitude to latitude on the sphere mapped to the icosahedron:
as described at https://en.wikipedia.org/wiki/Latitude#Authalic_latitude . This is how I understand what you mean by:
Which would result in the only area discrepency (as computed on the WGS84 ellipsoid) being the pentagonal zones having 5/6th the area of the hexagonal zones. Thanks again so much, your input is invaluable to getting this all right and achieving our goals of DGGS interoperability. |
@jerstlouis That’s all I was getting at. As long as you’ve densified the edges in ISEA space and verified the area is stable you should be fine.
This is a deep topic beyond our scope here. I suggest reading the discussion on using DGGS to represent DGGS cell geometries (including vertices) in Sahr 2019, in particular the section entitled “EXACT REPRESENTATION AND SUBSTRATE GRIDS”. Using cells to represent “vector” points can be applied to any vector use case, not just the representation of cell vertices. The last time I looked at the OGC standard it appeared to me that you were claiming to solve the problem of unifying raster/vector datasets by rasterizing the vector datasets, which is not a solution to the problem at all. But we are indeed off topic…
That looks to me like you’re doubling down on DGGS cell/zone indicating an areal unit only.
I don’t know enough about OGC terminology to know what your definition of “gridded coverage cell” is. Just taking the words as english I would say that DGGS units are “gridded coverage cells”.
Could be. I don’t think geodetic latitude is defined the same on a sphere and an ellipse, but I’m not a geodesist.
I’m not sure what point you’re trying to make. Why would we use the authalic radius otherwise? DGGRID does not currently support an ellipsoid model, but implementing the authalic conversion has been on my todo list since the beginning. But we’re getting off topic here….
I believe so, yes. |
|
Thanks again @sahrk .
That is one approach (which is what the DGGS-JSON representation does for example), which I would have thought offers some partial solution to the problem, if not an ideal / perfect solution. There are other approaches, possibly better suited for vector data, such as clipping vector geometry to a particular areal zone (and generalizing the geometry with a minimum distance treshold between points based on that zone's resolution), or expressing the geometry vertices as zone IDs representing their centroids. I'll have to read your 2019 paper -- thanks!
That's how the terms are specified in Topic 21 definitions, so I'm not really "doubling down", but at least for OGC standards we have to try to align with that terminology. I think addressing the use cases of representing / addressing "points" is a different topic, regardless of the terminology used.
I believe that OGC "gridded coverage" as defined by Topic 6 (e.g., see 5.3.2 Grid on page 11-12) implies a rectilinear grid, which DGGS units are not necessarily. So for gridded coverage cells, although the cells themselves are not necessarily shaped in a particular way, they are at least organized in a rectilinear grid. e.g., the rectilinear pixels/gridded coverage cells in this PNG zone data response: https://maps.gnosis.earth/ogcapi/collections/gebco/dggs/ISEA3H/zones/A6-0-E/data.png do not correspond to finer zones/cells of ISEA3H (this output is in EPSG:4326, but it would also apply to an ISEA projection output CRS). I can see however how the definition of "Grid" in Topic 6 could be extended to be more flexible so that most DGGs could potentially qualify. Whereas for a DGGS-JSON representation of the same zone data, individual values for sub-zones correspond to finer hexagonal zones (as visualized above).
I actually always wondered about that!!! :) Since it seems that we've been doing data integration from WGS84-referenced data in GNOSIS and PYXIS without using authalic latitudes. Now I understand the purpose of that authalic radius because of that authalic latitude concept, and I'll have a go at implementings this WGS84 geodetic latitude => authalic latitude conversion. |
|
@jerstlouis I fear we've hijacked this pull request but just wanted to add one more thing that might interest you...
This is the basic way to encode vector locations. This approach encodes the location at a single precision, and then (possibly) indexes that encoding using a hierarchical indexing. That is the best that is possible using systems like H3 or A3HT. But some hierarchical indexing systems allow for multi-precision quantization. See Sahr 2013 or the corresponding slide deck. |
|
Sorry for diverging into the above long parenthesis in this PR, but the discussions were highly insightful and did address the questions I had. To recap and summarize and try to move this forward:
Thanks! |
As far as I understand, what is not currently supported is the However, there are no tests for any of these even for the forward projection, which at least in TDD practice would imply they do not work :) I question whether Actually, there is one test expecting a failure: I imagine it's testing a resolution too high? |
me neither. So perhaps just mentions that the inverse direction is only available when none of those options is specified?
yes, this just tests we don't run into undefined behavior found by a fuzzer which generates "random" PROJ strings, and was fixed per 4c8a5cb |
This is actually confusing and probably wrong. The ISEA projection is strictly on the icosahedron. As discussed above, the paper confusingly talks about several different projections on different solids, but ISEA is specifically about the icosahedron. The original code had some hints that the author wanted to implement different solids, but only the icosahedron is implemented.
This does not read well at all to me, and the description says nothing about the discrete global grids that this projection is widely used for. I'm proposing to change the description to:
and add this note regarding the use of ellipsoids:
|
jerstlouis
force-pushed
the
master
branch
6 times, most recently
from
ab841ac
to
e66e3b0
Compare
- Integrating code ported from Java to eC, then eC to C++ implementing the inverse ISEA projection - Originally from Franz-Benjamin Mocnik's ISEA implementation found at https://github.com/mocnik-science/geogrid/blob/master/src/main/java/org/giscience/utils/geogrid/projections/ISEAProjection.java (MIT License) - NOTE: The inverse only supports the default planar options
- Eliminating several redundant trigonometric function calls - Single shared definition of dodecahedron vertices (centers of icosahedral triangular faces) for forward and inverse - Shared data structure for forward and inverse - Consistent use of 0-base index for triangular faces
|
What is this if (pj_param(P->ctx, P->params, "trescale").i) {
Q->radius = ISEA_SCALE;
}How do these |
- NOTE: There may have been an undocumented +rescale= parameter which will no longer work
cf https://github.com/OSGeo/PROJ/blob/master/src/param.cpp#L112 |
|
The That Those The isea_ptdd() is eerily familiar, looking very much like the 60 deg rotated / 30 deg horizontally sheared 5x6 ISEA space I was talking about above: convert projected triangle coords to quad xy coords, return quad number isea_rotate(pt, downtri ? 240.0 : 60.0);
if (downtri) {
pt->x += 0.5;
pt->y += cos30;
}which can be achieved with an affine transform: Does all this really belong in PROJ? It seems like all that code belongs in discrete global grid reference system libraries like geogrid and DGGRID, not in a projection library. With the exception of what This None of this is tested, and none of it is documented properly. Would there be any way to know if anyone is using these things, and if so, could we suggest to move this code to a separate module that can build on top of the default (Does PROJ have a deprecation mechanism?) |
Ask first on https://lists.osgeo.org/mailman/listinfo/proj , which I assume will be met with silence :-), and then remove them and wait for a few months to see if people complain when the new PROJ version hits them. In any case, that should be done in a separate pull request
We don't really have a runtime warning level. If we wanted to do a stage deprecation, perhaps the isea.rst page should be updated to mention that parameters X, Y, Z will be removed in PROJ 9.6 Besides that, is there anything left for this PR? It looks good to me as it |
Are you OK with the undocumented I did separate commits in case you wanted to pick and choose what to merge. Feel free to squash them if you prefer.
That would be a good idea. Would you want to do this separately and/or after asking on the mailing list?
Just let me fix this |
separate (clean) commits is fine
I would (let you) asking on the mailing list first, and in the absence of negative feedback, proceed to update the doc |
- Avoid calling isea_rotate() in isea_tri_plane() to rotate coordinates in bottom triangles by 180 degrees.
Thanks, will do. I just tried adding tests for |
I would appreciate if we could merge quickly this PR. Currently I need to click on a button in github UI each time you update the PR, so that github actions are run, because, as you don't have any commit merged yet into this repository, github doesn't trust you yet... |
Sorry to be spamming PR updates. |
|
I've just merged this PR. |
- Additional `azimuth` parameter for polyhedron orientation which is necessary to place the second vertex e.g., for R. Buckminster Fuller's Dymaxion Orientation ( see https://webpages.sou.edu/~sahrk/dgg/orientation/dymorient.html ). Latitude and Longitude of first vertex is not enough. - Clarify that to reference spatial information to the DGGRS, geodetic latitudes on the WGS84 ellipsoid must go through a geodetic latitude ==> authalic latitude conversion ( see https://en.wikipedia.org/wiki/Latitude#Authalic_latitude ) to map them to the authalic sphere. - See related discussions in OSGeo/PROJ#4211
+orient=poleor+orient=isea, but+orient=polewas not tested)